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Abstract. - We present an analysis of the statistical properties of hydrodynamic field fluc- 
tuations which reveal the existence of precursors to fingering processes. These precursors are 
found to exhibit power law distributions, and these power laws are shown to follow from spatial 
g-Gaussian structures which are solutions to the generalized non-linear diffusion equation. 



Over the past fifteen years, a nonextensive variant to the Boltzmann-Gibbs formulation 
of Statistical Mechanics has been developed for the analysis of systems away from the equi- 
librium state [1] characterized by the appearance of non-exponential distributions and power 
laws. For instance, for diffusion type processes, non-exponential distributions are obtained 
from generalized Fokker-Planck type equations, and it has been shown, in a generic manner 
from a generalization of classical linear response theory, that such distributions may arise 
from first-principle considerations [2] . The goal of this letter is (i) to show that in fingering 
phenomena, patterns are preceded by precursor processes, and (ii) to present an analysis of 
the statistical properties of hydrodynamic field fluctuations in these precursors; we find power 
law distributions which are shown to follow from the solution of the generalized non-linear 
diffusion equation. Although we do not use nonextensive statistical mechanics to explain our 
results, we find mathematical similarities which suggest a possible connection as discussed in 
the concluding paragraph. 

The phenomena investigated here arise before the onset of fingering, a generic phenomenon 
that results from the destabilization of the interface between two fluids with different mobilities 
in systems such as a shallow layer or a porous medium, when the fluid with highest mobility is 
forced through the medium filled with the other fluid. Here before means that the constrained 
fluid is in a state where no fingering pattern is as yet visible, but where hydrodynamic field 
fluctuations are enhanced as precursors to the onset of fingering. The analytical form of the 
statistical properties of these precursors are compatible with the solution of the generalized 
diffusion equation [2] which has formally the same structure as the "porous media equation" 
[3], but where the diffusion coefficient depends on the solution of the equation. This leads 
to the fact that the diffusion process is classical in the sense that there is linear scaling with 
time, but the solutions are not Gaussian: they have the canonical q-exponential form [1]. 

To date - at least to the best of our knowledge - such precursor properties have not yet 
been obtained from laboratory measurements. Here we use two methods: (1) a mesoscopic 
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Fig. 1 - Illustration of lattice Boltzmann (LBGK) simulation of viscous fingering (flow direction 
downwards). The mixing length between the two fluids is clearly identified by the shadow zone 
around the concentration = 0.5 contour line (color code indicating the concentration of the invading 
fluid from red (Cj = I) to blue (Ci = 0)). 



approach, the lattice Boltzmann simulation method [4], which is based on a kinetic theoretical 
analysis where the macroscopic description is not pre-established, and (2) a phenomenological 
approach, Darcy's law [5] which is solved numerically. Our results show that (i) the system's 
non-linear response to the driving force produces dynamical behavior in accord with the gen- 
eralized diffusion equation while no such phenomenological prescription was injected a priori, 
(ii) the statistics of the hydrodynamic fluctuations reveal the existence of precursor processes 
before fingering patterns become visible, (iii) the distribution of these field fluctuations are 
given by power laws (with an exponent \q\ < 1), and (iv) these power laws follow from space 
integration of two-dimensional g-Gaussian structures. 

These structures form a landscape of alternating upward and downward (/-Gaussian "blobs" 
(see Fig. 2 below) which are precursors to the fingers, and are reminiscent of vortices with 
alternating parities in 2-d turbulence For instance, experiments in turbulent Couette- 
Taylor flow [6] have shown that data obtained from quantities extracted from velocity field 
measurements exhibit g-exponential distributions which have been analyzed with theoretical 
arguments based on nonextensive statistical mechanics [7] . Other physical systems have been 
shown to exhibit similar type distributions interpreted as a consequence of nonextensivity [1] 
or superstatistics [8]. We suggest that in the early stage of the fingering process, the signature 
of nonextensivity can be found in the recasting of interacting Gaussian structures into a sum 
of independent q- Gaussian blobs. 

We consider a two miscible fluids system confined between two plates with narrow gap, a 
configuration known as the Hele-Shaw cell [9]. The Hele-Shaw geometry is intrinsically three- 
dimensional, but the effective destabilization of the flow can be described in the 2-d plane, 
given that the flow has a Poiscuillc profile in the third dimension. So Hele-Shaw flow can 
be simulated with the two-dimensional LBGK equation, the Bhatnagar-Gross-Krook single 



(^Here the Reynolds number has low value, but the relevant control parameter, the Peclet number, is high. 
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Fig. 2 - Concentration fluctuations field showing landscape of g-Gaussian "hills and wells" (flow 
direction South-East); simulation data with concentration contour lines, color code indicating highest 
positive values (red) to largest negative values (magenta) of c. 



relaxation version of the lattice Boltzmann equation [4]. One emulates 3-d flow by introducing 
a drag term, thereby simulating a system with a virtual cell gap in the third dimension: the 
drag enters the LBGK equation as a damping term. The simulation method has been described 
elsewhere [10]: essentially the system consists of a 2-d box with vertical length L y = 2048 
nodes, and horizontal width L x — 1024 nodes, with horizontal periodic boundary conditions; 
initially the box is filled with fluid 2 which is then displaced by fluid 1 injected uniformly from 
the top of the box. Fluid 1 (with damping coefficient (3\ < P2) invades the system through 
a constant force applied along the y-direction, which produces a flow with Peclet number 
Pe ~ 200. 

The initial concentration profile is a step function, which because of mutual diffusion, 
develops into an erfc y, on top of which white noise is imposed (along the x-direction) to 
trigger the instability. When the instability develops, it produces a moving fingering pattern 
as illustrated in Fig.l. Here we consider the early stage where fingers are not yet visible 
(they have not developed as in Fig.l shown here for the purpose of illustration), just below 
the dynamical transition where the exponent [i of the growth of the mixing length of the 
interfacial zone, L m i X oc i M , changes from /i = 1/2 (diffusive regime) to fi ~ 2.3 (see Fig. 3 
in [10]). In this early time regime, the flow produces local concentration gradients which 
induce mobility fluctuations thereby triggering vorticity fluctuations. 

We measured the following local fluctuations: 
• concentration field (fluid 1): c(x,y) = C(x,y) — C(x,y); 

■ transverse velocity field: v x {x,y) — V x (x,y) — V x (x,y); 

■ longitudinal velocity field: v y {x,y) — V y (x,y) — V y (x,y); 

■ and vorticity field: to{x,y) = d x v y (x,y) — d y v x (x 7 y), 

where the overline denotes average along the ^-direction. The concentration field in Fig. 2 
shows a "landscape" of alternating hills and wells. In each blob, the concentration field 
exhibits a two-dimensional g-Gaussian profile as illustrated in the left panel of Fig. 3 obtained 
by a section plane cut through the extrema in Fig. 2 parallel to the x-axis. The corresponding 
concentration distribution, P(c), is obtained by taking the values of c(x,y) measured over all 
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Fig. 3 - (a) Left panel: Concentration fluctuations profile (solid line) obtained from the simulation 
data by a section plane cut through the hills and wells in Fig.2; analytical function Eq. J2J (dashed 
line): the convex and concave curves are connected by summation of four g-Gaussians; (b) Right 
panel : Concentration fluctuations distribution (dots = data) with power law (solid line): P(c) — 
A\c\- q with q = 0.73. 



space, and is shown to follow a power law, P(c) oc \c\~ q with q = 0.73 (see right panel of Fig. 3). 
Note that the power law exponent stays the same over the entire range of the distribution. The 
other local fields exhibit similar behavior and the corresponding distributions P(v x ), P{v y ) and 
P(u>) all show power laws (with an exponent < 1) which follow from space integration (see 
below) . 

Processes stemming from the generalized diffusion equation [2] (or from the generalized 
Fokker-Planck equation provided some conditions are met for the drift and dissipative func- 
tions [11-13]) can produce solutions of the g-exponential type. In the present simulations, 
we observe that the concentration field fluctuations exhibit a q-Gaussian profile - while the 
system is still in the diffusive regime (~ t 1 / 2 ) - which suggests that the concentration c(r,t) 
is governed by the generalized e?-dimensional non- linear diffusion equation [2] ( 2 ) 

(1) 

A general solution to Eq.Q is indeed a q-exponential of the form [2] 

c(r) = co e-T« |r|2 = Co [1 - (1 - q) lq \r\ 2 ]^ , (2) 

where c$ and 7 g are time-dependent quantities, and q + a = 2. This solution (taken at fixed t) 
is used to perform a fit to the simulation data, and the concentration profile found from the 
simulation results has exactly the form of as shown in Fig. 3a. 

On the other hand, the velocity field is related to the concentration through Darcy's law 
(see e.g. [9]) which can be written for the stream function ip (in two dimensions) as 

(d d "N 1 

( 2 )On the r.h.s. of Q one can rewrite formally D d r c a as D a d r c, where D a = a D c a 



dc(r, t) 

Ft 



1 d 

r d— 1 Q r 



r d - x D{r)^c a {r,t) 
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Fig. 4 - (a) Left panel: Stream function distribution (dots = data from numerical solution of Eq.©) 
with power law (solid and dashed lines) P(tp) oc \^\~ q © with q' = 0.83 for \ip > 0|, and q' = 0.88 
for \ip < 0|. (b) Right panel: concentration field distribution (dots = data from numerical solution 
of classical advection-diffusion equation coupled to Eq.JHJ) with power law (solid and dashed lines) 
P(e) oc \c\' q with q = 0.78 for |c > 0|, and q = 0.82 for |c < 0|. 



with k x = dc/dx, and R = 91n/3/i9c, where (3 is the damping function which, in the fingering 
simulations, controls the drag and depends on the concentration c. Equation (0), coupled 
with the classical advection-diffusion equation for c, was solved numerically under conditions 
similar to those used for the lattice Boltzmann simulations. For both the concentration field 
and the stream function, wc obtain q-Gaussians in agreement with the data obtained from 
the LBGK simulation results; in particular, for the stream function, we have 

r/>{r) ^^ [l-(l-q')^r 2 }^ , (4) 

where ipo is a normalization constant 

We now show that the distribution which follows from the q-Gaussian spatial profile is 
a power law. Generalizing the ^-exponential in |0J to an arbitrary power A of the radial 
variable r, we have 

V(r) = e"^' rX = [1 - (1 - q')^}^ . (5) 

(for simplicity we have incorporated the normalization into ip), and the corresponding prob- 
ability distribution function P(ip) reads (in d dimensions) 

/>oo 

P(ijj)=w d r d - l \dr\5{^{r)~xP), (6) 
Jo 

where Wd is an integration constant. Then J5j is inverted to give r x = ^(l-g') = — 4T7 m ?' ^> 
which is used in ijfJJl to obtain 

PW = ~* = [- In, ^ r< • (7) 
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Fig. 5 - (a) Left panel: Vorticity field showing landscape of g-Gaussian "hills and wells" (flow direc- 
tion South-East); simulation data with vorticity contour lines, color code indicating highest positive 
values (red) to largest negative values (blue) of u). (b) Right panel: Contour lines of vorticity field 
(simulation data) superimposed on color concentration fluctuations field (simulation data); color code 
indicates highest positive values (red) to largest negative values (blue) of c. The inflection line of the 
concentration field (C ~ 0.5) is about two thirds down from top of picture (flow direction downwards). 



For two-dimensional systems (d = 2) with g-Gaussian spatial profile (A = 2), we obtain the 
power law 

Pty) = (8) 

2(p q > 

which is exactly what we find for the stream function distribution shown in Fig. 4a where 
a fit to the simulation data gives an exponent q < 1 (indicating statistical importance of 
large values of the stream). The distribution P(c) (Fig. 4b) is obtained along the same lines 
and yields the value q ~ 0.80 ± 0.05 to be compared with the value 0.73 ± 0.05 obtained 
from the lattice Boltzmann simulations (see Fig. 3). Singularities (elliptic singularities [14]) in 
the distribution function occur where the field has extrema (corresponding to the maximum 
of the absolute value of the blobs). According to Eq.0J, the value of q obtained from the 
concentration profile and distribution data yields a ~ 1.25 in the diffusion term of Eq.Q, 
thus an effective diffusion coefficient D a tx c 25 (see footnote (2)). These results indicate that 
the generalized diffusion equation governs effective transport in fingering phenomena. 

We also measured the vorticity field which, as shown in Fig. 5a, exhibits a pattern of al- 
ternating vortices which can be viewed as a landscape of hills and wells of radial g-Gaussians 
distributed along the concentration profile of the interfacial zone. Figure 5b shows the respec- 
tive locations of the vorticity field (contour lines) and concentration field (color coded) . 

Two questions arise from the present results: (i) what is the physical picture of the pre- 
cursors dynamics? and (ii) given that many systems exhibiting similar features have been 
analyzed with the nonextensive approach [15], we pose the reciprocal question: can one sug- 
gest a physical interpretation of nonextensivity from the structure of the precursors ? 
(i) We observe that, before any fingering pattern becomes visible, the driving force triggers 
pre-transitional fluctuations in the concentration field in the form of blobs organized spa- 
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tially according to the wavelength of the forthcoming fingers (Figs. 2 and 3a). Concentration 
fluctuations modify locally the C-dependent drag coefficient thereby inducing shear. As a 
consequence, the initially uniform velocity field is perturbed, and the resulting velocity fluc- 
tuations produce vortices with alternating parities and whose spatial sequence matches (with 
phase shift) the concentration fluctuation structures (Fig. 5). The configuration shown in 
Fig. 5b is reminiscent of the fingering vortices depicted in a schematic in [16] indicating that 
the mechanism described in [16] is already present in the precursors. 

(ii) The system studied here is representative of a generic class of driven noncquilibrium 
systems with spatial fluctuations, where q-exponentials and power law distributions are the 
signature of long-range interactions, and whose dynamical behavior is governed by advection- 
diffusion equations. What we have shown is that during the onset of fingering, one can 
identify precursors which exhibit features with mathematical properties typically encountered 
in nonextensive statistics. This similarity suggests a possible interpretation of the mechanism 
of nonextensivity. The driving force produces a spatial sequence of alternating structures, 
which, if they were independent, would exhibit ordinary Gaussian profile originating from 
local diffusion centers (5-functions), and would be described by classical diffusion. However, 
when growing, these Gaussians overlap, and what the analysis shows is that the structure 
formed by these overlapping Gaussians can be cast into a sum of scale invariant g-Gaussians 
with finite support (q < 1), i.e. these g-Gaussians are independent as corroborated by the 
fact that the power law distribution is shown (see Eq.©) to follow from space integration of 
a single ^-exponential. 

This work was supported by a grant from the European Space Agency and PRODEX 
(Belgium) under contract ESA/14556/00/NL/SFe(IC). 
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